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Transport in a Bacterial Porin 
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We have studied the dynamics of chloride and potassium ions in the interior of the OmpF porin 
under the influence of an external electric field. From the results of extensive all-atom molecular 
dynamics simulations of the system we computed several first passage time (FPT) quantities to 
characterize the dynamics of the ions in the interior of the channel. Such FPT quantities obtained 
from MD simulations demonstrate that it is not possible to describe the dynamics of chloride and 
potassium ions inside the whole channel with a single constant diffusion coefficient. However, we 
showed that a valid, statistically rigorous, description in terms of a constant diffusion coefficient D 
and an effective deterministic force F e g can be obtained after appropriate subdivison of the channel 
in different regions suggested by the X-ray structure. These results have important implications for 
popular simplified descriptions of channels based on the ID Poisson-Nernst-Planck (PNP) equations. 
Also, the effect of entropic barriers on the diffusion of the ions is identified and briefly discussed. 

PACS numbers: 05.60.Cd, 02.50.Fz, 87.16.Vy, 87.15.Vv 
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I. INTRODUCTION 

Understanding transport of charged solutes across the 
cell membrane is a problem of paramount importance 
in biophysics, since it is crucial to regulate many cell 
functions [l[ . This task is undertaken by transmembrane 
channels, so the characterization of the permeation of 
charged particles, ions in particular, is essential. In addi- 
tion, recently there has been a renewed interest to com- 
prehend transport through biological nanochannels for 
possible applications in biotechnology @. 

Bacterial porins are macromolecular proteins located 
at the outer membrane that enable the diffusion of small 
molecules through the lipid bilayer. As porins are well 
characterized, both structurally and functionally, they 
represent model systems to study transport through bi- 
ological nanochannels. In particular, the ionic transport 
properties of the OmpF porin have been extensively stud- 
ied since the determination of its X-ray structure [|| . 

The simplest description of the complex problem of 
ionic transport across channels is based on the classical 
diffusion-advection transport equations: 



dc(z,t) dJ(z,t) 



dt 



dz 







J(z,t) = -D 



( dc{z,t) F(z) 



\ dz 
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Here, J(z,t) is the flux of particles of a given ion of 
charge q and c(z,t) its concentration profile, D is a 
constant diffusion coefficient characterizing the brown- 
ian motion (with a value in principle different from bulk 



diffusion coefficient) and F(z) is a deterministic force. In 
the classical Nernst-Planck equations, the deterministic 
force arises from the mean-field electrostatic potential, 
F(z) = -qd<t>{z)/dz whereas in more advanced ap- 
proaches it may contain steric or entropic forces 

In some cases, integrated solutions of Eqs. (|1|2|) 
under suitable boundary conditions and appropriate 
assumptions are employed. The Goldman-Hodkin-Katz 
equation or the Henderson equation are well-known 
examples, being extensively used in the literature to 
understand the conductance, reversal potential, and 
other measurable quantities of nanochannels [l| . In other 
cases, it is solved self-consistently in combination with 
Poisson's equation of electrostatics, forming the so-called 
Poisson-Nernst-Planck (PNP) equations. This approach 
allows the analysis of electrically charged nanochannels 
in a consistent fashion 0-EI1- In all cases, a continuous 
mean-field approximation is assumed, i.e. ions are 
treated not as discrete entities but as continuous charged 
densities that represent the space-time average of the 
microscopic motion of the ions. The diffusion coefficient 
of ions, D, characterize their motion within the channel, 
representing a spatial average over the region considered. 
In analyzing experiments, diffusion coefficients are often 
given bulk values when the comparison with experiment 
is qualitative or are treated as free parameters to adjust 
by fitting experimental data 
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In recent years, the increase in computer power and the 
development of new algorithms have permitted the study 
of ion transport through certain nanochannels emplo yin g 
all-atom Molecular Dynamics (MD) simulations |13 - ll5j . 
Such analysis captures all the atomic detail of the system 
(including the water molecules present), which in some 
cases is necessary to properly account for the interactions 
of the ions with the interior of the channel. However, 
these MD simulations still remain computationally very 
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expensive, so its use is hindered by limitations in the size 
of the system (~ 10 — 100 nm) and the duration of the 
processes (~ 100 ns) that such calculations are capable 
of tackling. Nevertheless, even for the analysis of larger 
systems or longer processes, all-atom MD simulations can 
be used to obtain parameters that coarser (less detailed) 
models employ in their description of the system. In par- 
ticular, MD simulations can be used to obtain diffusion 
coefficients of ions in the interior of nanochanncls, which 
are central input parameters in many theoretical descrip- 
tions of ion permeation (like, e.g., Brownian Dynamics 
simulations [Ril - tla or approaches based on the Nernst- 
Planck equation [20} ). Tieleman and Berendsen 21 char- 
acterized the diffusive dynamics of water molecules in the 
OmpF porin by using the Einstein relation for the mean 
square displacement (MSD). They calculated local diffu- 
sion coefficients which they averaged out over slices per- 
pendicular to the axis of the pore to obtain the variation 
of the diffusion coefficient along the axis of the channel. 
A similar approach was used with the OmpF porin [22[ 
and other channels [H, [H[ to determine the dependence 
of ionic diffusion coefficients along their axial direction. 
Such approaches, however, do not provide global diffu- 
sion coefficients to describe the entire channel, and it is 
not clear how they can be determined. To circumvent 
that limitation, Hijkoop et al. 24|, inspired on the work 
of Munakata and Kaneko [25| , used a method based on 
a first passage time (FPT) analysis to compute a global 
diffusion coefficient of water in the interior of an artifi- 
cial (uncharged) OmpF porin. Their analysis, however, 
ignores the influence of free energy barriers and wells on 
the local movements of water molecules at a specific po- 
sition along the axis of the pore. 

In this paper, we use a method based on a first pas- 
sage time (FPT) analysis of the particles' trajectories 
to study the transport of potassium and chloride ions 
through the OmpF porin. This procedure allows us to 
obtain the effective diffusion coefficients of the ions in the 
channel from all-atom MD simulations [U, [H, H3- O n 
the one hand, from the trajectories of the ions provided 
by the simulations one can calculate several FPT quanti- 
ties that characterize their passage through the channel. 
Note that such quantities are produced taking into ac- 
count all the complexity of the all-atom description. On 
the other hand, it is possible to derive analytical expres- 
sions of such FPT quantities for a model based on the ID 
diffusion-advection equations with constant parameters, 
Eqs. (|1I2[) . The comparison of the FPT quantities ob- 
tained from the MD simulations with the analytical for- 
mulas give answer to two relevant questions for the char- 
acterization of the channel. First, it determines whether 
the one-dimensional description of the dynamics of the 
ions in the channel based on Eqs. (11 121) is valid. Second, 
in case such a description is accurate, we can obtain the 
diffusion coefficients and effective forces that best define 
the channel by fitting the FPT quantities calculated from 
the MD simulations with their corresponding analytical 
expressions. 



The rest of the paper is organized as follows. In sec- 
tion [IT] we describe the structure of the OmpF porin and 
discuss the model and the details of the MD simulations. 
In section HIIl we define the first passage quantities which 
will be considered and give the corresponding analytical 
expressions for the ID model based on Eq. (|1I2[) . In sec- 
tion |IV] we provide the results of the analysis. We discuss 
the FPT quantities obtained for the whole channel and 
also study the FPT properties of different regions within 
the channel. From this analysis, effective diffusion coeffi- 
cients and forces are obtained for the different regions. In 
section |Vj the conclusions of the study are given, as well 
as a table with the obtained diffusion coefficients that 
characterize the ionic transport of the channel. 



II. SYSTEM AND SIMULATION METHODS 

The biological nanochannel considered here, the outer 
membrane protein F (OmpF porin), is a transmembrane 
channel located at the outer membrane of the bacteria 
Escherichia coli. The structure of the WT OmpF channel 
is illustrated in Figure [T] The crystallographical struc- 
ture of the wild type OmpF channel is well known from 
long ago [1, Hll . It is made of three identical monomerical 
nanopores, each one assembling into a large 16-stranded 
antiparallel /3-barrel structure enclosing the transmem- 
brane pore (Figure [T]) . Each aqueous pore has a diam- 
eter between 1-4 nm, being constricted around half-way 
through the membrane by a long loop. It is a relatively 
wide channel in the sense that it allows the simultaneous 
permeation of both cations and anions (in hydrated form) 
and also the passive diffusion of larger dipolar molecules. 
The transport of monovalent ions and molecules across 
this channel is known to be dominated by electrostat- 
ics, as demonstrated by studies considering mutations of 
the channel (with different pore sizes and different elec- 
trostatic charges) and different permea ting spe cimens of 
different sizes and charge distributions |l9l l29l . |30| . 

In this paper we study the transport of potassium and 
chloride ions through the OmpF porin using a FPT ap- 
proach from results of all-atom molecular d yna mics sim- 

"13. For the 



ulations previously published in Rcfs. [151 . 
sake of completeness, we summarize here the most impor- 
tant features of the simulation. The system is composed 
of the OmpF porin inserted in a POPC membrane in 
contact with a 1M KC1 aqueous ionic solution, see FigfS] 
The whole system is thermalized at T r = 296K. The co- 
ordinates of the OmpF porin arc obtained from its X-ray 
structure ||. The protonation states employed for the 
titatrable residues of the protein are the same as those 
specified in Ref. [33j . which results in an overall charge 
of —lie per monomer. The molecular dynamics simula- 

We 
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tions were carried out with NAMD version 2.6 
used the force-field for protein-lipid simulations supplied 
by NAMD 2, which is a combination of the CHARMM22 
force field, parameterized to describe protein systems, 
and the CHARMM27 force field[35j], parameterized to 
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0.62 nm 



FIG. 1: (Color online) Protein structure derived from its X- 
ray structure (a) Top view of the OmpF porin. (b) Cross 
section of a single monomer, showing its dimensions. Figure 
produced using VMD [12] 



describe lipid systems in aqueous media. Water was also 
modeled after the CHARMM force field [IH, which de- 
scribes water using a slightly modified version of the stan- 
dard TIP3P model. Within this force field, ions are de- 
scribed as charged Lennard- Jones spheres. All simulation 
details can be found in Refs.fTol. l3l|. 

To study the ionic transport properties of the chan- 
nel, the whole system is under the action of a uniform 
electric field perpendicular to the lipid bilayer (along the 
negative z axis, sec Figj2]). The magnitude of the field is 
Eext — 14.22mV/nm, which creates a potential drop of 
~ 200 mV across the membranc+channel system. From 
a computational perspective, the high electrolyte con- 
centration, of approximately 1M, is convenient, since the 
large number of ions included in the simulation signifi- 
cantly helps to obtain statistically meaningful averages 
of the ion properties inside the nanochanncl. Such high 
concentrations, however, are not unrealistic. Similar con- 
centrations are often used in experiments to measure 
ion fluxes and selectivity properties of ion channels 
The total simulated time is ~ 30ns, considerably longer 
than previous MD simulations done of the OmpF porin 
[T3 . 122I ] , which also helps to improve the statistical aver- 
ages of ion properties inside the channel. 



III. FIRST PASSAGE TIME (FPT) APPROACH 

Let us consider a diffusing Brownian particle under 
the action of an external force in a one dimensional in- 
terval [0,L]. To characterize its dynamics in the interval 
one can compute first passage quantities, which deal with 
the event of the particle reaching for the first time any 
of the boundaries (first passage event). The central mag- 
nitude that characterizes these events (from which other 
first passage quantities can be calculated) is the distri- 
bution of first passage times, given a starting position of 
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FIG. 2: (Color online) Snapshot from the MD simulation, in 
which the dimensions of the simulation box have been added. 
In the picture, two slabs of 1M KC1 aqueous solution are 
separated by a POPC lipid bilayer in which the trimer OmpF 
porin is inserted. To favor the passage of ions, an electric 
field along the negative z axis is applied. For clarity, only a 
fraction of the lipids of the membrane are represented. Figure 
produced using VMD [32| . 



the particle within the interval. This is obtained from 
the probability that the particle reaches for the first time 
any of the boundaries at a given time t. This first pas- 
sage problem is, thus, approximately equivalent to the 
problem of a diffusing particle in an interval [0,L] with 
absorbing boundaries (2(| [27} . 



A. Calculation of FPT quantities from MD results 

The analysis based on a first passage approach is well 
suited to characterize the dynamics of the ions in the 
channel from the results of molecular dynamics simula- 
tions. One can easily calculate the different first-passage 
quantities characterizing the ions dynamics inside a given 
interval within the channel from inspection of the indi- 
vidual trajectories provided by the simulations' results. 

A first passage quantity that we employ in our anal- 
ysis is the Hitting Probability. It is defined as the prob- 
ability that a particle which starts at an initial posi- 
tion zq reaches for the first time either the lower bound- 
ary at z = 0, P-(za), or the upper one at z = L, 
P+(z ) = 1 — P_(z ). The dependence of P±{zq) on z 
is determined by the presence of free energy barriers and 
external forces acting on the ions. The calculation of the 
hitting probability from the MD trajectories of the ions is 
simple. At each timestep of the simulation we identify the 
ions contained within the layer [z$ — A/2, zq + A/2], with 
A = 0.1 nm (thincr layers were tested with identical re- 
sults). We then trace all those particles until they reach 
for the first time any of the virtual boundaries of the 
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interval, moment at which they loose their identity (al- 
though they are not physically removed from the system) . 
By keeping track of the number of ions that reach each 
boundary we can compute, once we have gone through 
all the simulation, the hitting probabilities of each type 
of ion for a given initial position zq. The same procedure 
is repeated for different starting points zq (separated by 
A) until the whole interval is covered to obtain the hitting 

probabilities from MD simulations, P± D '(zq). 

Another magnitude of interest in the application of the 
FPT approach to the analysis of ionic dynamics in the 
channel is the Mean Exit Time (MET), T(z ). It is de- 
fined as the average time taken for a particle to reach 
for the first time any of the boundaries given that its 
initial position is Zq. To calculate the MET from the 
MD simulations results we divide the system into layers 
of width A as done before, and at each timestep identify 
the ions within the layer centered at zq. Those parti- 
cles are tracked until they reach either boundary and the 
time elapsed is recorded. From the collection of residence 
times of the particles starting at zq we compute the av- 
erage time to reach either boundary. This is repeated 
for all zq in the interval to obtain the MET from MD 
simulations, T( MI) )(z ). 

The other FPT quantity that we use in the charac- 
terization of ion dynamics in the channel is the Sur- 
vival Probability, S p (t). It gives the probability at time 
t of finding a particle within the domain [0,L] that will 
cross the domain from one boundary to the other with- 
out crossing the entrance at intermediate times. For the 
calculation of this quantity, at each timestep of the sim- 
ulation we identify those ions that enter the considered 
interval and, for each particle, we record the time of en- 
trance into the interval and a tag which specifies the 
boundary of entrance. As done in the previous cases, we 
track those particles until they reach any of the bound- 
aries, where they loose their identity. If the boundary 
of exit is the same as the boundary of entrance the tra- 
jectory is discarded. If the particle exits through the 
opposite boundary, the exit time is recorded so we can 
compute the time taken by that particle to cross the in- 
terval. From the collection of the crossing times it is then 
straightforward to calculate the survival probability from 
MD simulations, S^ ID) (t). In the presence of an exter- 
nal effective force dragging the ions, there are many more 
events of ions crossing the interval in the direction of the 
dragging force than in the opposite direction. In order to 
provide statistically meaningful results, only the survival 
probability corresponding to the passage of ions in the 
direction of the dragging force will be considered. 



late the corresponding theoretical FPT quantities 
P2 h) (z ;{a n }),T( ih )(z ;{a n }),Si th \t;{a n }), where 
{a n } are the parameters of the model. By comparing 
such theoretical quantities with those computed from 
MD simulations, P^ D) (z ),T^ MD )(z ),S^ MD \t), we 
can test the applicability of the model and also, under 
favorable conditions, determine the parameters {a n } 
that define such model. Note that the parameters 
thus obtained contain in an effective fashion all the 
information on the geometry and interactions of the 
all-atomic simulation. 

Here, we consider the simplest diffusion model based 
on ions characterized by a constant diffusion coefficient D 
moving under the influence of a constant effective force 
F e g. For such a simple model, the classical transport 
Eqs. (|II2p apply and simple analytical expressions for 
the FPT quantities can be deduced in terms of the pa- 
rameters D and F e g. These analytical expressions are 
derived by noting that the one-dimensional FPT problem 
is approximately equivalent to the problem of a diffusing 
particle in a ID interval with absorbing boundaries. For 
detailed derivations, the reader is refereed to standard 
textbooks such as Ref. [27|. Here we limit ourselves to 
list the results relevant for our analysis. 

The Hitting Probabilities are given, in terms of D and 
F eS , by g3 



P { . h \z ) 



Note that P± h \zo) depend solely on one of the param- 
eters of the model, F e e. Consequently, the comparison 
of Eqs. © with the results from the MD simulations, 

P±* D \zo), is well suited to obtain the effective force 
acting on the ions. 

Within the considered theoretical model, the Mean 
Exit Time reads (H 



-F B f t z /k B T 



-F ctt L/k B T 



1 _ e -F cit L/k B T 
~l _ g — -Faff z /k B T 
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(3) 
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Finally, the analytical expression for the survival prob- 
ability, Sp h (t), can also be deduced from the diffusion- 
ad vection equations, Eqs. (|1I2[) . [27j : 



(-1) 



n-l„2 



-a„t 



^ n 2 + (F eS L/2irk B T) 2 



(5) 



B. Analytical results for a simplified model 

The FPT quantities obtained from the MD sim- 
ulations can be used to check the validity and 
to extract the parameters of a certain theoretical 
model. Based on such a model, one can calcu- 



with a n = (D/L 2 ) [n 2 n 2 + (F cB L/2k B T) 2 ]. At times 

t 3> l/o!2, Sp h \t) is dominated by the first term, with 
cti = (D/L 2 ) [ir 2 + (F cS L/2k B T) 2 ], which directly de- 
pends on the diffusion coefficient. 

In the analyses done below, based on the comparison of 
the FPT results from the MD simulations with Eqs. ([3]), 
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(O and ([4]), the interval length L needs to be considered 
slightly longer than its actual length in the MD simu- 
lations [13]. This extra distance, which appears to be 
dependent on the type of particle and the interval, has 
been studied and identified before as the Milne extrap- 
olation length [ID, The necessity for such a (very 
small) correction emerges because the absorbing bound- 
ary conditions assumed in [27| to derive Eqs. ©, (jU) and 
(|5j) are not exact for the virtual interval considered in the 
analysis of the simulations. 



IV. RESULTS AND DISCUSSION 

A. Test case: Bulk Electrolyte 

To check the validity of the method described above, 
we analyzed first the dynamics of the ions in bulk, 
far from the channel and membrane, which we know 
is describable using a ID Nernst-Planck equation 
fEqs. (|H2[ )') with constant parameters and whose results 
can also be obtained from the analysis of the mean 
square displacement (MSD). 

For this analysis we have selected a 2nm-thick slab of 
electrolyte far enough from the channel and membrane, 
— 6.0nm< z < —4.0 ran (see Fig. [5]). The delimiting 
planes of the slab perpendicular to the z-axis act as vir- 
tual absorbing boundaries, while in the transverse direc- 
tions (x, y) periodic boundary conditions are applied. 

The protocol that was followed to obtain the param- 
eters of the theoretical model, D and F e g, is as fol- 
lows. The effective force parameter F e s is determined 
first from fitting the hitting probabilities of reaching the 
upper boundary of the specified region, pf MD >( Zo ),with 
Eqs. The dependence of these probabilities on the 
initial position of the particle is determined only by the 
effective force F c g, being a straight line for F c ff = and 
increasing its curvature as the parameter grows (concave 
for F c ff > and convex for F e g < 0). Second, the dif- 
fusion coefficients are obtained from the comparison of 
the MD results for the survival probability, Sp MD \t), 

with Eq.© at t > l/a 2 . The slope of log[S^ MD \t)} di- 
rectly gives the diffusion coefficient, provided that F e g is 
known. An alternative way to obtain the diffusion coeffi- 
cients that is used to check the consistency of the results 
consists of the comparison of the computed mean exit 
times, T^ MD ^ (zo), with the analytical expression Eq. (QJ. 

From the calculation of the hitting probability 
p(MD) we nave verified that, due to the much higher 
resistance of the mcmbrane+channel, in the bulk there is 
a negligible drop of electrostatic potential. The effective 
force F e ft obtained from the comparison of P| M (zo) 
with Eq. ([3]) , which is determined by the value of the 
electric field, is not distinguishable from (see Fig. 

The analysis of the survival probability provides the 
diffusion coefficients for the CI - and K + ions. From the 
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FIG. 3: Hitting probability, P+(zo), of reaching the upper 
boundary of the bulk region (located at z = — 4nm) for (a) 
chloride ions (b) potassium ions as a function of their ini- 
tial position, 2o . The linearity of the curves is an evidence 
supporting that F c a — 0. 



slope of the fitting straight line represented in Figs. Q we 
obtain Dqi = 2.11nm 2 /ns and Dk = 2.05nm 2 /ns with 
the help of Eq. ©. 

These values compare very well with the values ob- 
tained from the study of the mean exit time. By fitting 
the results for the mean exit time computed from the MD 
simulations, 7 1 ( A/D ) (zo), with its corresponding analyti- 
cal expression, Eq. we obtain Dei = 2.08nm 2 /ns and 
Dk = 2.07nm 2 /ns. Such agreement of the parameters, 
obtained from two independent methods, shows the con- 
sistency and robustness of the FPT procedure employed 
to extract the values of the diffusion coefficient D and ef- 
fective force F e g that characterize the simple ID model. 
An extra independent test was done to check the reliabil- 
ity of the results obtained from a first passage time ap- 
proach. Simulations of a 1M KC1 electrolyte in a box with 
periodic boundary conditions were performed and calcu- 
lated the diffusion coefficients using Einstein's formula 
for the MSD [3?|]. The resulting diffusion coefficients are 
Dqi = 2.07nm 2 /ns and Dk = 2.08nm 2 /ns, which are in 
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FIG. 4: Survival probability in bulk of 
potassium ions. 
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FIG. 5: Mean exit time in bulk of (a) chloride ions, (b) potas- 
sium ions as a function of their initial position Zq. 
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very good agreement with both sets of values obtained 
from a first passage time approach. The results are also 
in good agreement with the experimental bulk diffusion 
coefficients for K + and Cl~ in a 1M KC1 concentration 
at 298K: D^ p) = 2 nm 2 /ns and D<£ xp) = 1.85 nm 2 /ns 

m. 



B. Channel 

In the following, we present the results for the First 
Passage quantities corresponding to the transport of ions 
across the whole channel. The channel is defined as the 
central part of the OmpF trimer in which the monomers 
have a rigid beta barrel structure, leaving out the more 
flexible loops and turns located at the edges of the pro- 
tein. In our simulation, the channel approximately ex- 
tends from z = —1.4 nm to z = 2.1 nm, see Fig. 
The results for the FPT quantities, Figs. (|6|) and (O, are 
obtained from the analysis of the trajectories of the ions 
of the all-atom MD simulation of the system, which con- 
tains the exact geometry and considers all interactions 



within the CHARMM force field. They provide valu- 
able information on the ion transport inside the channel. 
From Fig. ([6]) we can infer the existence of an energy bar- 
rier located in the vicinity of z ~ 0.5 nm that prevents 
the passage of both ions, since the probability of crossing 
z ~ 0.5 nm is very small. The results on the MET rep- 
resented in Fig. ([7]) provide the characteristic residence 
times for both ions as a function of their position within 
the channel. From those values we can obtain the average 
time of residence in the channel for both ions, Tci = 822 
ps for chloride and T k = 935 ps for potassium ions. The 
different position of the maximum of the MET for both 
ions, located at z < zl/ 2 for chloride ions and at z > zl/2 
for potassium ions (with zl/ 2 being the middle of the 
channel), reflects the overall action of the external elec- 
tric field (directed towards negative z-direction), which 
drags cations along the field and anions in the opposite 
direction. 

Another interesting outcome of the analysis is that the 
results for the hitting probabilities and mean exit times 
shown in Figs. © and ([7]), contrary to the bulk case, 
cannot be satisfactorily fitted by the analytical expres- 
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FIG. 6: Hitting probability, P+(zn), of ions reaching the 
upper boundary of the channel (located at z = 2.1nm) as a 
function of the initial position of the particle, zq. 




sions Eqs. and (jl]). Our results show that it is not 
possible to assign constant diffusion coefficients to the 
ions inside the channel. Therefore, this result cast some 
doubt on the adequacy of one-dimensional models based 
on the Nernst-Planck equation with constant diffusion co- 
efficients to describe the ion transport through the whole 
OmpF porin. 
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In order to reach a better understanding of the ionic 
transport through the OmpF porin, wc divided the sys- 
tem into different regions and, through the comparison 
of the FPT quantities, P^ MD) (z ),T ( * MD '> (z ), S { p MD) (t), 
with the corresponding analytical expressions, deter- 
mined the parameters that characterize the ion dynamics 
(i.e., the effective diffusion coefficient D and the effective 
force F ff) in each region. Based on its structure, we 
split the channel into three different regions: Constric- 
tion, defined as the central region between z = —0.3 nm 



FIG. 8: (Color online) Different regions of the system: 
Vestibulel, —1.4 nm< z < —0.3 nm. Constriction, —0.3 
nm< z < 1.1 nm. Vestibule2, 1.1 nm< z < 2.1 nm. Fig- 
ure produced using VMD [321 ]. 



and z = 1.1 nm; Vestibulel, defined by z = —1.4 nm and 
z = —0.3 nm; and Vestibule2, in between z = 1.1 nm and 
z = 2.1 nm (see Fig. ©). 

C. Channel Constriction 

We employed the same protocol used in the analysis of 
the ions in bulk to study the dynamics of K + and CI" in 
the constriction region. Although the quantitative agree- 
ment between the first passage quantities computed from 
the simulations results and the corresponding analytical 
expressions is not accurate, the general functional form 
of the MD results is well captured by the analytical ex- 
pressions with the appropriate constant parameters (see 
Figs. ©, (|TU)) . ([TTjO . Such an agreement ensures the ad- 
equacy of the ID diffusion-advection, Eqs. (UJ [2]), to de- 
scribe the dynamics of the ions inside the constriction. It 
must be stressed here the reduced number of ions present 
in the constriction, which results in poor statistics and 
might contribute to the observed deviations. 

From the analysis of the hitting probabilities, we ob- 
tain the values of the effective force acting on the ions (see 



Fig. ©). The effective force is F ( 



-1.9 k^Tr/nm 

and F^S = 1.6 kgT r /nm for potassium and chlo- 
ride ions, respectively. As expected, potassium ions are 
dragged along the direction of the external electric field, 
whereas chloride ions are dragged in the opposite direc- 
tion. The slightly different magnitude of both forces is 
consistent with the observation that anions and cations 
follow different pathways in the constriction, as has been 
reported in previous MD studies 22j and in experiment 
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|39| . Assuming that the section of the channel along the 
z-axis is uniform in the constriction (see Fig. ©), it 
follows (see below) that the effective force acting on the 
ions is determined solely by the local electric field act- 
ing on them. From the values of the effective forces we 
can estimate the magnitude of the local electric field for 
potassium ions, \F^ ] /e\ = 0.049 V/ nm, and for chloride 
ions, \Fj£ l) /e\ = 0.041 V/nm. 
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FIG. 9: Hitting probability, P+(zn), of ions reaching the up- 
per boundary of the constriction (located at z = l.lnm) as a 
function of the initial position of the particle, zq. The concav- 
ity of the curve in the chloride ion case indicates the action 



of a force along the z-direction, F e 



(co 



1.6 k s T r /nm. The 



convexity of the curve in the potassium ion case indicates 



the action of a force in the opposite direction, 
ksT r /nm. 
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FIG. 10: Survival probability of ions in the constriction. 



1200 

1000 

& 800 
E 

t 600 

X 

HI 

c 

2 400 
2 

200 



CI": MD simulations • 

CI": Fit 

K + : MD simulations A 

K': Fit 

A A 




0.25 0.5 
z (nm) 



The fit of the time dependence of the survival proba- 
bilities obtained from the MD simulations Sp (t) with 
Eq. (|5|) provides the diffusion coefficients Dei = 0.18±0.1 
nm 2 /ns, D K = 0.21 ± 0.1 nm 2 /ns. As seen in Fig. (flQ]), 
the statistics in this case is poor, so the results for the 
diffusion coefficients are associated with large errors. 

These values compare well with the values obtained 
from the analysis of the mean exit time: Dei = 0.26 
nm 2 /ns, Dk = 0.22 nm 2 /ns. Also, as seen in Fig. (fTTT) . 
the profile of the mean exit time is well reproduced by 
Eq. (fj| with such parameters. Note that the diffusion 
coefficients of both ions in the constriction are roughly an 
order of magnitude smaller than the diffusion coefficients 
in the bulk. 



D. Channel Vestibules 

A similar procedure was employed to analyze the dy- 
namics of the ions at the vestibules of the constric- 
tion, defined as the regions between z — — 0.3nm and 
z = — 1.4nm (Vestibulcl) and z = l.lnm and z = 2.1nm 
(Vestibule2) (see Fig. flg)). 



FIG. 11: Mean exit time of ions in the constriction as a func- 
tion of their initial position. 



The first general observation is the remarkable good 
agreement between the first passage quantities computed 
from the MD simulations results and the corresponding 
analytical expressions with the appropriate constant pa- 
rameters (see Figs. (|T2"j) , (fT4|) and (p~5|) V This indicates 
that the dynamics of the ions in the vestibules is prop- 
erly described by the simple ID model based on Eqs. ()l|2p 
with the appropriate parameters D and i^ff- 

For the vestibules, the results for the hitting prob- 
abilities (zq) give the surprising results shown 
in Figs. (|T2"j ). In the two regions, the effective forcces 
acting on both anions and cations have the same sign, 
being negative in Vestibulel {F^ 1 ^ = — 1.4ksT r /nm for 
chloride ions, and F e ( ^ = -2.6k B T r / nm for potassium 
ions) and positive in Vestibule2 (-Fg ff = 5.8kBT r /nm 
for chloride ions, and F^g = 4.0kg T r /nm for potassium 
ions). This indicates that the dominant contribution 
to the effective force dragging the ions is not charge 
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FIG. 12: (a) Hitting probability, P+(zo), of ions reaching the 
upper boundary of Vestibulel (located at z = — 0.3nm) as a 
function of the initial position of the particle, zo- For chloride 



ions, 

7 (K) 



^eff ' = ~~ 1.4ksT, /nm is used. For potassium ions, 
F^g' = -2.6k s T r /nm. (b) Hitting probability, P+(z ), of 
ions reaching the upper boundary of Vestibule2 (located at 
z = 2.1nm) as a function of the initial position of the particle, 
zo- For chloride ions, F^ 1 ^ = 5.8ksT r /nm. For potassium 
ions, F^P = 4.0k B T r /nm. 



dependent and its direction must be along the z-axis 
(and of larger magnitude) in Vcstibule2 and in the 
opposite direction in Vestibulel. Also, the magnitudes 
of the effective forces acting on each ion are different. 
The effective force acting on potassium ions is greater 
than the force on chloride ions in Vestibulel, but it is 
smaller in Vestibule2. 

Such results correspond well with the action of charge 
independent entropic forces, which arise when the dy- 
namics of particles through channels is described with 
one dimensional models. In treating channels of varying 
cross section, the analysis of the dynamics of ions can be 
reduced, under appropriate plausible hypothesis, from a 
three-dimensional problem with a complex boundary to 



a one-dimensional problem of diffusion along the channel 
axis d-d, H3, E3 ■ Assuming that the distribution of par- 
ticles in the planes perpendicular to the axis of the chan- 
nel is close to its equilibrium distribution, one obtains the 
modified Ficks-Jacobs equation. Following Reguera and 
Rubi's formulation Q, in this approximation the flux of 
particles, Jfj, can be expressed as: 



J FJ {z,t) = -D 



dc(z,t) 1 dA(z) 



()-: 



kuT r dz 



c(z,t) 



(6) 



Here c(z, t) is the ID particle concentration (with dimen- 
sions of particles per unit length), A(z) = U(z) — TS(z) 
is a free energy given by the electrostatic energy U(z) = 
—qEz and the entropy S(z) = fcfllii/i(z), where h(z) 
is the dimensionless cross section of the channel. This 
entropic term accounts for the variation of the space ac- 
cessible for the diffusion of the particles, dragging them 
towards positions z of maximum entropy. Note also that 
Eq.© is a particular case of Eq. (|1I2[) with a deterministic 
force arising from free energy barriers, F = —dA(z)/dz. 

From Eq. ^ we can infer that the effective forces ob- 
tained in the vestibules from the analysis of first passage 
quantities can be decomposed into two additive contri- 
butions: 



F e fi = Fei + F e , 



(7) 



Here, F e i is the component of the force arising from elec- 
trostatic interactions, whose sign depends on the charge 
of the ion under consideration. F entr accounts for the ef- 
fective entropic force, which depends solely on the shape 
of the channel and which points towards the direction 
that maximizes the available space. 

With the help of Eqs. ([7]), we can obtain the two con- 
tributions, F e i and F entr , from the values of the total 
effective force computed from the first passage analysis, 



F, 



(K/Cl) 



F, 



off 



F, 



(Cl) 

ir 



(8) 



In Vestibulel, the resultant effective electrostatic force 



is F. 
F, 



(K/Cl) 



trostatic component is F^^ Cl ^ 



^0.6 ksT r /nm, and the entropic part is 
2ksT r /nm. In Vestibule2, the resultant elec- 

^0.9 ksT r /nm, and 
the entropic part is F entr = 4.9 kgT r /nm. 

To confirm the existence of a charge independent en- 
tropic force arising from the variation of the available 
space along the axis of the channel, we have also ana- 
lyzed the first passage properties of water molecules at 
both vestibules (see Figs. (fl~3|) ). In this case, the only con- 
tribution to the effective force is entropic. For Vestibulel 
we find F entr = —1.5 kgT r /nm, and for Vestibule2 we 
find F en tr = 4 ksT r /nm. 

At this point we can interpret qualitatively the results 
obtained from the analysis of the hitting probabilities at 
the two vestibules. In Vestibulel, the effective section 
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FIG. 14: Survival probability of ions in (a) Vestibulel (b) 
Vestibule2. 



of the channel decreases along the z-axis, which induces 
the appearance of an entropic force in the same direction 
as the external electric field (antiparallel to the z-axis). 
Thus, for cations both contributions to the effective force 
have the same direction whereas for anions they have op- 
posite directions. As the entropic part is dominant, the 
effective force on both ions is directed towards decreasing 
z values, although its magnitude is larger on potassium 
ions than on chloride ions. In Vestibule2, on the contrary, 
the effective section of the channel increases along the z- 
axis, so the entropic force is directed towards increasing 
z values, and the two contributions of the effective force 
have the same direction for anions and opposite direc- 
tions for cations. Hence, in this case the total effective 
force is in the direction of the z-axis and its magnitude 
is larger for chloride ions than for potassium ions. 

The entropic force can be independently estimated us- 
ing Eqs. © from the dependence of the effective radius of 
the channel on the z-coordinate, as computed from the 
MD simulations. We obtain an average F entr = —1.1 
ksT r /nm at Vestibulel and an average F entr = 3.1 
ksT r /nm at Vestibule2. The discrepancies in the pre- 
cise values obtained for the entropic forces were ex- 



pected, since the expressions used are only valid for 
smoother changes of the section of the channel along the 
z-coordinate. 



Once the effective forces have been determined we 
can obtain the values of the diffusion coefficients. The 
logarithmic fit of the time dependence of the survival 
probabilities, Sp MD \t), provides the diffusion coeffi- 
cients Dei = 0.65 nm 2 /ns and Dk = 0.66 nm 2 /ns in 
Vestibulel, and D C i = 0.49 nm 2 /ns, D K = 0.33 nm 2 /ns 
in Vestibule2 (see Figs. (|14[ 0. Independently, from the 
analysis of the mean exit times, we obtain Dei = 0.57 
nm 2 /ns, Dk = 0.59 nm 2 /ns in Vestibulel and Dei = 
0.44 nm 2 /ns, D K = 0.31 nm 2 /ns in Vestibule2 (see 
Figs. (US])). Again, such an agreement of the parameters, 
obtained from two independent methods, shows the con- 
sistency and robustness of the FPT procedure employed 
to extract the values of the diffusion coefficient D. 
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FIG. 15: Mean exit time of ions in (a) Vestibulel (b) 
Vestibule2 as a function of their initial position. 



V. CONCLUSIONS 

In this paper we have studied the diffusive dynamics of 
chloride and potassium ions in the interior of the OmpF 
porin under the influence of an external electric field di- 
rected along its axis. From the results of extensive all- 
atom molecular dynamics simulations of the system we 
computed several first passage time quantities to char- 
acterize the dynamics of the ions in the interior of the 
channel. Such quantities, calculated considering all the 
complexity of the system, provide valuable information 
on the residence times and presence of free energy barri- 
ers for ions in the interior of the OmpF porin. These FPT 
quantities were also used to test the validity of simplified 
ID descriptions of the dynamics of the ions in the inte- 
rior of the channel based on the use of constant diffusion 
coefficients for the ions. 

Overall, our results show a slow and complex dynamics 
for the ions inside the channel. The sluggish ionic dynam- 
ics inside the channel is characterized by large residence 
times (~ 0.8 ns for Cl~ and ~ 0.9 ns for K + ) and very 
large mean exit times with peaks around ~ 2 ns, see Fig- 



TABLE I: Diffusion coefficients in the different regions. Two 
values are provided, one inferred from the analysis of the sur- 
vival probability and the other from the study of the mean 
exit time. 





Diffusion coefficient D (nm^/ns) 


Region 


Ion 


Survival Prob. 


Mean Exit Time 


Bulk 


K+ 

cr 


2.05 
2.11 


2.07 
2.08 


Constriction 


K+ 

cr 


0.21 
0.18 


0.22 
0.26 


Vestibulel 


K+ 

cr 


0.66 
0.65 


0.59 
0.57 


Vestibule2 


K+ 

cr 


0.33 
0.49 


0.31 
0.44 



urc[71 which indicate a strong interaction of the channel 
with the ions. The complexity of the dynamics is demon- 
strated by the fact that it is not possible to characterize 
chloride and potassium ions in the interior of the channel 
with a single diffusion coefficient. However, it is possible 
to find well-defined diffusion coefficients for some regions 
of the protein channel. 

Based on its X-ray structure, we have divided the chan- 
nel into three different regions: constriction, Vestibule 1 
and Vestibule 2. Employing the same FPT analysis we 
demonstrate that each of such regions is describable by 
a simple model characterized by a constant diffusion co- 
efficient D and a constant deterministic force F e g acting 
on the ions. The diffusion coefficients D which charac- 
terize the dynamics of chloride and potassium ions in the 
different regions, obtained by two independent methods, 
are given in Table |U Let us remark that these diffusion 
coefficients vary by almost an order of magnitude from 
the bulk to the interior of the constriction. The physical 
origin of the effective force F e g is also different in the 
different regions of the channel. The effective forces on 
the ions at the constriction, as inferred from the FPT re- 
sults of the MD simulations, correspond to the presence 
of an intense local electric field directed along the applied 
electric field which drags potassium ions along and chlo- 
ride ions in the opposite direction. In the vestibules, the 
effective forces observed on the ions are consistent with 
the expected effect of entropic forces which arise from 
the variation of the available space along the axis of the 
channel (see Eq.©). 

We believe that our results have important implica- 
tions for the modeling of ionic channels. In this field, 
there is a great interest in having simple, ID models, 
with the capability of predicting the essential features 
of these extremely complex systems. Such models (like, 
e.g., the PNP approach) are based on the classical ID 
transport equations [l[ (see Eq. (|ll2[) ') and heavily rely on 
the existence of a well-defined diffusion coefficient. Our 
work strongly suggests that in general a direct use of 
a constant diffusion coefficient for the whole channel is 
not justified. However, it also shows the plausibility of a 
multi scaling approach, in which these models would play 
an important role. The idea is to take into account the 
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complexity of the structure of the protein by obtaining 
diffusion coefficients from detailed MD simulations in re- 
gions of different geometrical properties. The coefficients 
thus obtained can be used in the classical ID descriptions 
of ionic channels, such as the PNP approach. Also, we 
recall that our work further demonstrates that in such 
ID descriptions, entropic forces such as those described 
in Refs.JJd, HO, El[ must be taken into account. Work 
in this direction of a multi-scaling approach is now under 
way. 
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